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SUMMARY 


This project (NAG 9-574) was meant to be a three- year research project. However, due to 
NASA’s reorganizations during 1992, the project was funded only for one year. Accord- 
ingly, every effort was made to make the present final report as if the project was meant to 
be for one- year duration. Originally, during the first year we were planning to accomplish 
the following: (1) we were to start with a three-dimensional flexible manipulator beams 
with articulated joints and with a linear control-based controllers applied at the joints; 
(2) using this simple example, we were to design the software systems requirements for 
real-time processing, introduce the streamlining of various computational algorithms, per- 
form the necessary reorganization of the partitioned simulation procedures, and assess the 
potential speed-up realization of the solution process by parallel computations. 

The three reports included as part of the final report address: (1) the streamlining of 
various computational algorithms; (2) the necessary reorganization of the partitioned sim- 
ulation procedures, in particular the observer models; and (3) an initial attempt of recon- 
figuring the flexible space structures. We wish to state that much of the real-time effort 
via parallel computations will continue under a NSF grant as part of High Performance 
Parallel Computing Initiative « Other two aspects also constitute important attributes for 
real-time simulation of space operations. As such, they will also be pursued through other 
grants in the near future. 

We thank Dr. John Sunkel for encouraging us through the difficult period of what turned 
out to be a short, intense yet productive endeavor. 
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9.1 INTRODUCTION 

This survey is a follow-up on earlier ones (Park [1]; Park and Felippa [2]; Felippa and 
Park [3]; Park [4]) on direct time integration methods. The algorithmic characterization 
of the integration formulae offered therein, namely, stability, accuracy and implement- 
ation aspects, remains largely intact. Readers wishing to familiarize the algorithmic 
characterization may refer to the references dted above plus Hughes and Belytschko 
[ 51 . 

What we are about to cover herein reflects a steady shift of research thrusts in 
computational dynamics since the mid-1980s, from discipline-oriented dynamics to 
system-oriented dynamics, from sequential computations to parallel computations, 
and from efficiency/accuracy concern to system model improvements/refine- 
ments. The specific topics we survey in this chapter thus reflect their maturing stages; 
these developments do not fit into a coherent theory or categorization at the present 
time. 

Since time integration algorithms have been presented within the context of linear 
structural dynamics for most instances, first we report on computational methods for 
non-linear multibody dynamics. Major challenges in the development of computational 
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methods for multibody dynamics analysis have been the conservation of both energy 

and momentum, system constraint violations, and simulation speed. We will address 
some of these issues. 

The second topic we will report is methods for the solution of coupled-field 
problems, primarily methods for control-structure interaction (CSI) problems. The 
esign, modelling, analysis and real-time operation of CSI systems are one of the 
most intensely researched activities in recent years with applications ranging from 
aeroelastic tailoring, vibration control of reflectors deployed in low-earth orbits to 
active vibration control of suspension systems. The third topic we will present 'is a 
computational method for transient thermal-structure interaction problems This 
technique is relevant to the analysis of the thermal response of high-speed transport 
plane as well as integrated electronic chip thermal management problem. 


9.2 SOLUTION TECHNIQUES FOR MULTIBODY DYNAMICS 

The equations of motion for multibodies are characterized by two key features: 
hi? n ° n Inear kinematical relations and complex constraints. It is not the purpose 
o his chapter to make an exhaustive survey of available solution techniques. Rather 
we will examine selected techniques that meet our needs: computer implementability, 
a ap a ion to large-scale simulation, robustness and efficiency, in that order, 

ere are three aspects of solution techniques for multibody dynamics (MBD) 
an ysis. First, we must have at hand an efficient and accurate algorithm for updating 
the kinematical quantities such as angular orientations, angular velocities. Second, 
direct time integration of the equations of motion that correspond to the unconstrained 
states of multibodies must be performed. Third, an accurate and efficient treatment 
ot constraints is essential if the numerical solution is to maintain the given holonomic 
and nonholonomic constraints. In practice, the three aspects are intertwined so that 
one must achieve a careful balance in the employment of three computational aspects 
As computer implementation of the three require different strategies, we will discuss 
than separately first and bnng them together in the solution procedures. 

he numerical solution procedure for MBD systems which we describe herein is 
termed a staggered MBD solution procedure that solves the generalized coordinates in 
modul< ; from that for the constraint force (Park and Chiou [6]; Park et al. 
17, 81). A major advantage of such a partitioned solution procedure is that additional 
analysts capabilities such as active controller and design optimization modules can 
e easily interfaced without embedding them into a monolithic program. The solution 
o the equations of motion for constrained multibody systems, unlike typical structural 
dynamics problems, must satisfy at each time step the system constraints, whether 
holonomic or non-holonomic or time-specified manoeuvres. Because of this distinctive 

V He and COSt of a multibod y analysis package can be strongly 

\ ow efficiently and accurately the constraints are preserved during the 
numerical solution stage. S 

The system constraint forces can be either eliminated or retained depending upon 
the complexities associated with the elimination process. In general, it is preferable 
to eiumnare the constraint degrees of freedom if they are associated only with open 
ngic lirucs. On the other hand, if external torque or active control devices are 
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attached to those joints, it is computationally more advantageous to solve the 
constraint forces (or Lagrange multipliers) simultaneously together with the general- 
ized coordinates. Unfortunately, a straightforward way of computing the Lagrange 
multipliers can often lead to an unacceptable level of errors. The task of minimizing 
the propagation error due to violations of the constraint conditions is known as 
stabilization. We will describe a particular constraint stabilization which recasts the 
algebraic constraint conditions to a set of parabolic differential equations such that 
the constraint forces can also be integrated in time. 

To solve for the generalized coordinates of the multibody system, the equations 
of motion are partitioned according to the translational and the rotational coordinates. 
This sets the stage for an efficient treatment of the rotational motions via the 
singularity-free Euler parameters. The translational part of the equations of motion 
is integrated via a standard central difference algorithm. The rotational part is treated 
by a modified central difference algorithm in order to preserve the discrete angular 
momentum. Once the angular velocities are obtained, the angular orientations are 
updated via the mid-point implicit formula employing the Euler parameters. 

When the two algorithms, namely, the modified central difference algorithm for 
the rotational coordinates and the implicit staggered procedure for the constraint 
Lagrange multipliers, are brought together in a staggered manner, they constitute a 
staggered explicit-implicit procedure as detailed below. 


9 . 2. 1 Equations of motion for muttibody systems 

To motivate ourselves for the development of solution procedures for the multibody 
dynamics problems, let us introduce the following equations of motion: 


Md = Q - B r A, 



(9.1) 


<D(d, d) = 0 


(9.2) 


where M is the system mass matrix, d is the generalized velocity vector, u is the 
translational degrees of freedom, 0) is the angular velocity vector, B = d<b/d& is the 
constraint projection matrix, A is the constraint force vector, O are the constraint 
conditions that are imposed either on the subsystem boundaries or on the kinematical 
relations among the generalized coordinates, t is the time, ( ' ) denotes time 
differentiation, and Q is the generalized applied force plus non-linear inertia forces. 

We observe from (9.1) and (9.2) that the task for solving the governing 
multibody dynamical equations constitutes three computational procedures: accurate 
computations of the constraint force. A, or their elimination from the equations of 
motion (9.1), updates of the angular orientations, and the direct integration of the 
translational displacement u. To this end, we will first examine two distinctive ways 
of handling A. 
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9.2.2 Techniques for handling constraint conditions 

In principle, it is better to eliminate the constraint conditions, if possible, if the 
corresponding forces are not needed for design or interface with other analysis 
modules. For example, if the system consists of open-tree configurations and no 
active controller is applied, then it is best to eliminate the joint constraint attributes. 
On the other hand, when the system includes multiple closed-loop configurations or 
active controllers are present on several joints, then it becomes important to compute 
the Lagrange multipliers as accurately as possible. 

First, one can easily eliminate the system constraint forces via a coordinate 
partitioning strategy whenever any or all of the system components possess an 
open-tree topology. In the second procedure, we present a stabilization procedure 
for solving the Lagrange multipliers. A distinct feature of this stabilization procedure 
is that it can be implemented in a stand-alone module, thus can be interfaced not 
only with the equation solver for rigid-body systems but with that for flexible-body 
systems as well. 


Parallel Implementation of coordinate partitioning technique 

In this technique, a projection matrix that spans the null space of the constraint 
Jacobian matrix <D M is first constructed (see, e.g., Wehage and Haug [9]). A parallel 
methodology (Chiou [10]; Chiou el al. [11]) based on an arrowhead algorithm then 
can be applied to the resulting complementary set of equations of motion. We will 
present the procedure for open-tree systems. For a system that contains closed loops, 
a cut-joint technique can be used so that the present scheme can be equally applied. 

Let us introduce a projection matrix A such that, when its transposed matrix acts 
on the constraint force B T k, it gives 

A T B T A = 0. (9.3) 

This projection matrix can be obtained by expressing the total generalized velocity 

“ ~ <u T co T > in terms of the independent velocities d f and their time derivatives 
as 


d = AcF, d = Ad^ + Ad f . (9.4) 

Due to the property of (9.3), premultiplying the equations of motion (9.1) by A T 
yields 


A'Md = A T Q ( 9 .5) 

In conventional procedure, d in the above equation is replaced by (9.4b) and ~d! 
is then solved from the reduced equations of motion. In the solution to be described 
below, instead of solving the reduced equations of motion, we augment (9.4b) to 
(9.5) to form an arrowhead matrix equation: 


M -MAljdj _ fMAd'l 
-A T M 0 Jld'j j-A T Qj 


(9.6) 
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which can be partitioned as 
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(9.7) 


where n is the total number of bodies in the system. Decomposed in a manner 
convenient for parallel computations, one obtains 


M ; d, + D„. „ + n d' = g ; , ; = 1 n 

Ld,„ + ,, d ; = g' 

where 


(9.8) 


n n 

X ^(« + i./> = ~ X Ar Mr' «+ 1) = “ MrA' 1 = * ' •' 71 

;=■ i /=i 

n 

g, = (MAd') ; , j = 1 ", g' = - £ A T Q r 

/-t 

Each diagonal submatrix M ; represents the local mass matrix which is decoupled 
and can be factorized concurrently. An off-diagonal submatrix D (>ll+1) denotes the 
coupling between connecting bodies in the system. Since M is a constant matrix, 
(9.8) becomes 

d, = M-' (D„.„ +1) d' - & ). (9.9) 

Substituting (9.9) into (9.8b) gives a form of Schur complement : 

XD (n+ ,,M- , D () .„ +l) d^= Z D ( „ +w) M- l Q ; - ZD (n+1 >d'. (9.10) 

;=i ;= 1 ]= I 

The preceding treatment of the reduced equations of motion provides several 
parallel computational features. First, the parallelism can be exploited by mapping 
each processor onto a group of bodies so that independent computations such as 
the left-hand side of (9.10) can be performed concurrently. Second, since M is a 
constant mass matrix, it needs to be factored only once. Third, to solve for d f , a 
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parallel sparse solver may be utilized. Finally, once d' is obtained, computations of 
d from (9.4) is also easily parallelizable. 


Stabilization of constraint violations 

When the Lagrange multipliers cannot be eliminated or are to be retained for other 
purposes, one must solve for them. It has been known for some time that a 
straightforward direct time integration of the governing differential equation (9.1), 
augmented with a linearized form of the system constraints (9.2), often incur 
unacceptably high errors in the numerical solution. Of several techniques proposed 
to date, perhaps the method proposed by Baumgarte [12, 13] is the earliest known 
stabilization technique for computing the constraint forces. 

While the method due to Baumgarte works relatively well, it requires an a priori 
determination of stabilization parameters and the method breaks down when the 
number of independent system constraints change due to varying configurations. To 
cope with the varying system constraints without experiencing singularities, a 
penalty-based stabilization has been developed in Park and Chiou [6], The penalty 
procedure recasts the constraint equation in the form 


X = 



(9.11) 


as the basic constraint equations instead of (9.2) for both the holonomic and non- 
noionomic constraint conditions. We then time-differentiate, for the holonomic case, 
bhe above penalty-based equation to obtain: 



The numerical solution to the above companion differential equation is effected 
as roilows. J 

The constrained equation of motion (9.1) is integrated once using the implicit 
integration rule 


d— 1/1 = d" + <5d" +1/2 , 5=- 

l 


to yield 

d'~ L - = 5M-' , (Q" +,/I - B t A" +,/ 2 ) + d". (9. 13) 

This expression is substituted into (9.12) to obtain the stabilized differential equation 
tor rrse Lagrange multipliers 
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eX n + l/1 + c5BM _l B t X n + in = <5BM _I Q" + 1/2 + Bd". 

When the above equation is integrated once more with the trapezoidal rule, we 
obtain the following discrete update for X; 

(el + 5 2 BM- ! B t U" + i/2 = eX n + r^ l/2 (9.14) 

r^ I/2 = 5 2 BM-'Q n + ul + 5 Bd\ (9.15) 

The solution procedure for X presented above can now be invoked in a staggered 
manner in conjuction with the generalized coordinate solver to be described below. 


9.2.3 Time Integration of MBD equations of motion 

Once X is computed by the procedure in Section 9.22 or d when using the 
partitioning algorithm, one still needs to compute d, w and the angular orientations 
and their parameters at each time step. This task is carried out by employing an 
explicit-implicit transient analysis algorithm that exploits the special kinematical 
relationships of the generalized rotational coordinates vs. the angular velocity, 
namely, the Euler parameters. First, the integration of the translational coordinates 
and the angular velocity is accomplished by the central difference formula. It should 
be mentioned that the use of the central difference formula does impose a step-size 
restriction due to its stability limit (cu max h < 2) where co^ is the highest angular 
velocity of the system components for rigid-body systems or the highest frequency 
of the entire flexible members for flexible-body systems. The simplicity of its 
programming effort and robustness of its solution results can often become compelling 
enough to adopt an explicit formula, which is the view taken here. 


Explicit translational coordinate integrator 

In the conventional structural dynamics analysis, explicit time integration of the 
equations of motion by the central difference formula involves the following two 
updates per step: 


u" +,/2 = ii n - U2 + hu n 


-rr+l 


= u" + h u" +1/2 . 


(9.16) 


Unfortunately, the same integrator is not directly applicable to the rotational part 
of the equations of motion since (O is not directly integrable to yield the total 
rotational quantities except for some special kinematic configurations. This motivates 
us to partition d into the translational velocity vector, u, which is directly integrable 
and the angular velocity vector, a), which is not, and treat them differently, viz.: 
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( 9 . 17 ) 


The equations of motion can be written according to the above partitioning as 

"m, 0 
_ 0 



( 9 . 18 ) 


where 


Iq. 

Iq.' 


) fu — D u (u) — S M (u, q) — B 
K — DJ&) - SJu, q) — B 


iU \ 


( 9 . 19 ) 


in which the subscripts (u, m) refer to the translational and the rotational motions, 
respectively, f is the external force vector, D is the generalized damping force 
including the centrifugal force, S is the internal force vector including member 
flexibility, q is the angular orientation parameters, B u and B w are the partition of the 
combined gradient matrices of the constraint conditions (9.2). 

First, assume that u and q n+ 1/2 are already computed so that we can compute 
u and d> n+1/2 : 


{ u" +I/2 ' 
d >" +1/2 


= M 


-JQu 


Second, update the translational velocity at the step {n + 1) by 
u t+i = u» + ha n+l/1 . 

Third, we update the translational displacement, u, by 

u-3/z = u n+m + h{l n+i 


( 9 . 20 ) 


( 9 . 21 ) 


( 9 . 22 ) 


The updating of the angular orientations must be treated with care, and is described 
below. 


Updating of angular velocity via discrete angular momentum conservation 

hi order to update the angular velocity and angular orientations, we combine 
judiciously a momentum-conserving form of the central difference algorithm and the 
mid-pocnt implicit rule for computing the Euler parameters as follows (Park and 
Chiou [34D- 

First, we retain the basic central difference formula for computing the angular 
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velocity at the half steps: 

<o n+I/2 = (o n - xn + ha>" (923) 

where ft) is computed from the equations of motion, utilizing the angular velocity 
obtained from the discrete angular conservation law, as described shortly. 

Second, co n+ul is used to integrate the Euler parameters by 


q 


1 

2 


0 —ft) 
ft) —ft) 


q = A (m) q 


(9.24) 


where q — (<fo are the Euler parameters expressed in the body-fixed frame 

and d) denotes the skew-symmetric angular velocity tensor given by 


0 -ft) 3 co 2 

ft) = COj 0 —CO . 

— co 2 co 2 0 


(9.25) 


Implicit integration (9.24) by the mid-point rule yields 
q" +I - q" = hq n+1/1 = h A(ft)" +1/2 jq’ ,+1/i 
= -A(0)" +,/2 )(q’’ + I + q") 


(9.26) 


where A(o) n+,/:! ) can be viewed as the tangent matrix at the mid-configuration 
whereas q" +I/2 = (q n+1 + q")/2 is the mid-point average value. It was shown in Park 
el al. [7] that q n+I can be expressed as 


q" +I = 


I + - A(o)" +1/2 ) 


(9.27) 


where A = 1 + /i J (ft) 2 + co 2 + co*)/4. 

The updated Euler parameters q n+1 are then normalized to satisfy the constraint 
condition 


Once q" +I is available, the angular orientation at t" +1 is then obtained by 
R" +! = (2^ - 1)1 + 2qq T - 2^ 0 q, q = (, j 1 q 2 q J ) T 
where the superscript ( n + 1) is omitted on q. 


(9.29) 
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Third, the angular velocity <o" +1 that is needed to compute 6) for the next step 
integration is obtained via a discrete version of the angular momentum conservation 
law: 


= Jit "* 1 ' 2 ( 9 . 30 ) 

where M ra is the moment of inertia, (o is the angular velocity vector, and r is the 
applied moment, all expressed in body-fixed frame at the configuration t k ,0 <,k <n. 
For computational simplicity, we will choose k = 0, i.e. the initial configuration so 
that (9.30) becomes 


(R" +1 ) M w ar +I - (R") T M m a>" = h(R n + ul ) T r +vl 


(9.31) 


where the matrix R is the rotation transformation matrix from the inertial basis e to 
the body fixed configuration b according to 


b" = R" e 


(9.32) 


and co and t are now expressed in the superscript-indexed discrete b-bases. 
Therefore, the discrete angular momentum equation (9.30) becomes 


<o n 


= MJ'R" +1 | 


(R") M ra co" + /i(R" +,/2 ) T Q 


■) 


(9.33) 


where t is replaced by Q m from the Iefthand side of (9.18). 

It is noted that, whereas the standard difference formula (9.16) satisfies the linear 
momentum conservation for constant and linearly varying Q ^ (9.33) indicates that 
the use of a common basis is essential for the conservation of angular momentum. 
A similar approach was successfully utilized by Simo and Wong [15] in their 
development of a family of implicit algorithms. 

Fourth, the angular acceleration needed for the next time step (n + 1) is then 
computed tor each rigid body by 


d>'~ 1 =M-‘(Q" +i -tb" +1 M c 


,OJ 


(9.34) 


Equations 9.23)— (9.34) constitute the present modified central difference algorithm 
for integrating the rotational equations of motion in multibody systems. However, 
as many engineering multibody systems involve both holonomic and nonholonomic 
constraints, the computation of Q” + 1/2 is not as straightforward as (9.18) implies. 

For a momentum-conserving implicit algorithm, the reader may consult Simo and 
Wong [15’. For applications of the preceding MBD procedures to flexible multibody 
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dynamics, one may refer to Downer [16), Downer el al [17] and Downer and Park 
[18], who solve the flexible appendage deployment problem, among others. Finally, 
other MBD recent approaches can be found in Haug and Deyo [19], 


9.3 ALGORITHMS FOR CONTROL-STRUCTURE INTERACTION 
SIMULATION 

A second topic we should like to report in this survey is computational methods for 
the simulation of dynamic response of structures that are subject to active control 
forces. A general case of structural response under active control forces involves 
both large-angle rigid motions as well as transient flexible vibrations. Engineering 
examples include the manoeuvring of robotic arms, satellite attitude changes, 
deployment and vibration control of large space structures, and active vibration 
suppression of rotating machinery and vehicle suspension systems. 

When relatively small size models are adequate for describing the predominant 
motions and vibrations, the resulting active control strategies also can consist of a 
small number of actuators and sensors. However, as the structural model needs to 
be large due to the physical nature of the problem or due to the high-predsion 
requirement, so must be the size of the actuator/sensor numbers. It is for such large- 
scale control-structure problems that the following simulation methodology has been 
developed. 

Spedfically, simulation tasks for control-structure interaction (CSI) problems involve 
several computational elements and disdpline-oriented models such as structural 
dynamics, control law synthesis, state estimation, actuator and sensor dynamics, 
thermal analysis, liquid sloshing and swirling, environmental disturbances, and 
manoeuvring thrusts and torques. Because each of these computational elements can 
be large, it is usually not practical to assemble these computational elements into a 
single set of equations of motion and perform the analysis in its totality, which will 
be referred to as the simultaneous solution approach. First, the equation size of the 
total system can be simply too large for many existing computers. Second, the 
simultaneous solution by treating the coupled interaction equations as one system 
may destroy the sparsity of the attendant matrices, thus requiring excessive 
computations and storage space. Most important of all, any changes in the model 
or in the computational procedures will engender significant modifications of the 
required analysis software modules and hence require a painstaking software 
verification effort. 

In order to alleviate the aforementioned difficulties that exist in the simultaneous 
solution approach, a partitioned solution procedure that takes the following consider- 
ations into account has been developed. First, software development of any new 
capability is costly and time-consuming; thus, if at all possible, it is preferable to 
utilize existing single-field analysis modules to conduct the coupled-field interaction 
analysis. Second, the tasks for model generation and methods development of each 
field are best accomplished by relying on the experts of each single-field discipline. 
In order to accommodate both the software considerations and the single-field 
expertise, a partitioned (or divide-and-conquer) analysis procedure has been developed 
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r interachon anaJ ysis for direct output feedback systems (Belvin 
the rcT ^1 6 Vin P roce ^ ure abandons the conventional way of treating 

*d?Ke P co„S ”, r e "‘ ity nS,eai il ,rea,S lhe *«*"<' observed 

' COntroller/o J bserver interaction terms as separate entities. Thus, the CSI 

adoDted T r ^ 5 8*! 12ed 1 3S 3 coup ! ed 'fi e ^ P r °blem and a divide-and-conquer strategy 
adopted for the development of a real-time computational procedure. It should £ 

■ IOned tbat a similar concept has been successfully applied to other interaction 

s^s ems etal [22]l ^b-sfmctural interaction 

systems (Park [23]; Fehppa and Park [3]; Park and Felippa [2]), earth dam and oore- 

uid interactions (Park [23); Zienkiewicz et al. [24] and multibody systems ? with 
constraints (Park el al. [7]; Chiou [10]; Downer [16]). 

Firs^I^r pr0Cedure hin S L es on software and computational aspects. 

h ”? Cremen L t the e ^ uations ° f motion for each discipline are 

appHed^ S2 V by , C ° nS ' derin 8 the mteracKon terms as external disturbLres or 
! „ forC f' Second when necessary, computational stabilizaHon and accuracy 
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H 1 n ° u e that such P^iboned solutions of each discipline equaHon can 

be earned out on either a sequential or a parallel machine if certain message passS 
and memory-conflict issues are handled appropriately 8 P 8 


9.3. 1 Equations of motion for control-structure Interaction systems 

SntTby ST"’ ° f m0li ° n f ° r C ° nlro| - staKtoa ™y b. 


Structure: 

(a) 

M q "f f^q "h Kq — f + Bu + Gw 

qf°) = q<v q(0) = q 0 

Sensor output: 

(b) 

z = Hx + v 

Estimator: 

(c) 

x = Ax + Ef + Bu-|-L 7 
x(0) = 0 

Control force: 

(d) 

u = — Fx 

Estimation error: 

(e) 

y = z - (H,q + h„4) 


where 



and 
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A = 


0 

-M“ T K 


I 

— M~ l D ' 



FJ- 


In the preceding equations, M is the mass matrix, D is the damping matrix, K is 
the stiffness matrix, f(t) is the applied force, B is the actuator location matrix, G is 
the disturbance location matrix, q is the generalized displacement vector, w is a 
disturbance vector and the superscript dot denotes time differentiation. In (9.35b), z 
is the measured sensor output. The matrix H d is the matrix of displacement sensor 
locations and is the matrix of velocity sensor locations. The vector v is 
measurement noise. The state estimator in (9.35c) is assumed to be based on either 
the Kalman filter (Kalman and Bucy [26]) or a Luenberger observer [27] if the system 
is deterministic. The superscript ' denotes the estimated^states. The actuator output, 
u, is a function of the state estimator variables, q and q, and F x and F 2 are control 
gains determined for example by pole-zero placement or from the solution of an 
optimal control problem. The observer is governed by L, the filter gain matrix. For 
the special case where Lj is the null matrix (i.e. q = q), a second-order state estimator 
can be expressed as 

Mq + Dq + Kq = f + Bu + ML 2 y. (9.36) 

The effect of the above simplification on the observer stability and convergence is 
discussed in detail in Belvin [20] and Belvin and Park [28]. 


9 . 3.2 Simultaneous solution approach 

The numerical solution of (9.35) by the simultaneous solution approach begins with 
appropriate initial conditions, the feedback gain F and the filter gain L. The structure 
equation is written in first-order form 

x = Ax + Ef + Bu 4* Gw (9.37) 

where 


G = 


0 

M-’G 


The control gains and observer gains can be synthesized independently by noting 
that the stability of the structural system and the observer error stability are 
uncoupled. Introducing the error equation by the deterministic form (9.35) as 



e — x — x = 


(9.38) 
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and eliminating u yields 



(9.39) 


The stability of (9.39) is governed by the stability of [A — BF] and [A — LH1. 
Thus, the control gain F is suitably chosen from the matrix [A — BF] and the 
observer gain L from the matrix [A — LH], 

Subsequently, the simultaneous solution approach eliminates u and z from (9.35a c) 
and then solves the observer based closed-loop equations 



A -BF 

LH A-BF-LH 


{?} + {pif + i G l 

W (Ej (oj 


w. 


(9.40) 


The embedding effects of both the controller and the state observer result in an 
unsymmetric and non-sparse system matrix of dimension (4 N by 4N), where N is 
the number of structural degrees of freedom. Solution of (9.40) would require 
considerable software modifications of existing structural dynamics analysis programs 
tor large-scale CSI simulation purposes. In addition to losing the computational 
advantages associated with the finite element based CSI equation, the simultaneous 
solution approach requires the control law to be embedded into the observer model 
It the control law includes actuator, sensor and/or controller dynamics, additional' 
states must be added to the observer. This greatly complicates the observer model 
mu. . Sl 8 m ^ cant: software development for each class of control law dynamics, 

the difficulties associated with the simultaneous solution approach have prompted 
development of a partitioned solution approach for the CSI equations as described 


9.3.3 Stabilization for computations of control force and estimation 
error 

The partitioned solution procedure numerically integrates the structural equations of 
motion (9.3a a) and the observer equation (9.35c) by treating the control force u and 
e estimation error y as if they were applied terms in the right-hand sides. In this 

Wa L!r aHOn L ° f control - structure interaction systems using the partitioned solution 
procedure can be carried out by a judicious employment of three software modules: 
the structural analyser to obtain q, the state estimator to obtain q, and the stabilized 
solver for the control force u and the state estimation error y. Thus the partitioned 
procedure becomes computationally efficient and can preserve software modularity 
by exploiting the symmetric matrix form on the left-hand sides of (9.35a) and (9.35c). 

However, computations of the control force u and the state estimation error y by 
(9J.'a) ana 935e), respectively, can not only lead to an accumulation of errors but 
otter, can give nse to numerical instability. Hence, in order to make the partitioned 
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solution procedure robust, it is imperative to stabilize the partitioned solution process 

and/or numerically filter the solution errors in computing u and y. This is addressed 
below. 


First, we time-differentiate (9.35c) to obtain 
“ = F iq - F 2 q. 

Substituting q from (9.36) into the above equation, one obtains 
u + F 2 M-'Bu= - F 2 (M~’p + L z y) - F,q 
where the generalized rate of momentum p is given by 
p = (f - Dq - Kq). 


(9.41) 


(9.42) 


(9.43) 


The parabolic stabilization that led to equation (9.42) for computing the control 
law is sometimes called an equation augmentation procedure as it has not altered 
any part of the basic governing equation (9.35) except one time-differentiation of u 
assuming u exists. However, this assumption is later removed through time 
discretization as will be shown later in the chapter. 

, jVv ‘jSfjJlff the homo 8 en °us part of (9.42) has the filtering effect of the form 
[“ V* M in parlance of classical control theory, where s is the Laplace 

transform operator thus achieving the required stabilization. From the computational 
viewpoint, although F 2 M 'B is in general a full matrix, its size is relatively small 
as the size of u is proportional to the number of actuators placed on the structure’ 
bimilany, for the observer estimation error y one can stabilize its computation 
hrst by time-differentiating it 


7 + H„L 2 y = z - (H,p + H„q). 


(9.44) 


and substituting the observer equation into the above to obtain an augmented form 
ot the observer error equation: 


yH 0 L 2 y = z - H^l-'fp + Bu) - H rf q. 


(9.45) 


9 . 3.4 Stabilized partitioned equations and solution process 


The adoption of the second-order observer and the preceding stabilization thus 
replaces (9J5c , (9.35d) and (9.35e) by (9.36), (9.43) and (9.45), respectively, as 
summarized below. 
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Structure: 

(a) 

Mq + Dq + Kq = f + Bu + Gw 

q(°) = q 0 - q(o) = q 0 

Sensor output: 

(b) 

z = Hx + v 

Estimator: 

(0 

Mq + Dq + Kq = f + Bu + ML 2 y 

q(0) = 0, q(0) = o 

Control force.* 

(d) 

u + F 2 M _i Bu = - F 2 (M _, p + L 2 y) - F, 

Estimation error: 

(e) 

y + H„L 2 y = z - H^M -1 (p + Bu) - H rf q 


(9.46) 


6 J lffer wn e ^ etween t , he ori S inal governing equation set (9.35) and the 

the cf ? 3 \ r 9 ' 46 ‘ S ^ ° bstacle to com P u tot>'on of the control forces and 
trie state estimation error vector. 


9.3.5 Stability and accuracy ot partitioned solution procedure 

^ Ut f°Klf 1StabiIity , ana,ySiS ° f P artiHoned Procedures for a general coupled 
system .s still in an evolving stage. Hence, the analysis herein applies the relevant 

P^nTc Park r d B t in I2I]) in the P resent stabi ^ analysis ofX 
parthoned CSI solution procedure. The partitioned CSI solution procedure presented 

forinulS when discrehzed by unconditionally stable implicit time integration 
to obtain u^^ ald'^^T com Pff Hon fJ instability as it involves extrapolations 
procedure (nr The r ^ \ j' stablllt y analysis of the partitioned solution 

Et tn 5 P f dynamiCS ' ° bserver controller equations is 

l h j observer dlaraderisKcs H. L and th, controller 
ch^actenshcs B, F are specified Hence, the analysis that follows is restricted to an 

e aHon 5 ^ 61 "' y ~ ° J n wbat foll ° ws ' * * assumed that all of the stabilized 
equation set is time-discretized by a mid-point version of the trapezoidal rule. 

n or er o assess the computational stability of the present partitioned solution 
procedure, we construct a model single degree-of-freedom interaction equation as 

8 5tmCh,ral dampi " 8 3 m0dal Stn,rt “ al ^Hon of motion 


y + a) 2 y = - u 


(9.47) 


where y is a generalized coordinate and co is its associated frequency. 

tbe ™ ode ’ controller is assumed to consist of both the position and 
velocity feedback with appropriate weights given by F 

u = r\oi) y + C(O c y , 0J min < «, < ro ml , (9 48) 

?the Slirht 6 w dback which ranges from the minimum to the maximum 

JJv? s^uctural frequency contents, and rj and C are positive scalar coefficients that 
signify the strength of the position and the velocity feedback, respectively 

equ^n T 8 ^ ^ ** Stahilized form ° f (948) we have ^ model interaction 
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y + a) 2 y = — u 

u + (co^ = rjco 2 y - Cco c (o 2 y. ^ 9 ' 49 ^ 

Thus, the model interaction equations given by (9.49) represent the case of full state 
feedback. They do not, however, reflect the mode-to-mode coupling that can occur 
in reduced-order feedback controller. Nevertheless, an analysis of the computational 
stability using the above model interaction equations should shed insight on the 
overall stability of the present partitioned solution procedure. 

Time integration of the above model problem (9.49) by the mid-point rule 

= x" + Sx n+1 ' 2 

= i" + Sx n+U1 (9.50) 

_ 2x' ,+,/z — x" 

with y = 0 yields 


x n+ 1/2 

x" +ul 


f p +ul = ^ + 5yr 

(1 -I- <5{a>>; +1/2 = (t]0) 2 c — 5 C(o c a> 2 )y" p +m + (a> c y" 

(1 + 5 2 co 2 y +in = - 6 V- 1 ' 2 + y" + <5y" (9 ’ 51) 

y" +I = 2y' , + 1/2 — y", y" +in = (/ +1/2 - y”)/<5, y” +I = 2y" +I/2 - y" 
u n+1 = t]co 2 y +l + 


where y ” p * 1/2 is a stable predictor that is needed to initiate the staggered solution. 

Computational stability of the above difference equation can be assessed by 
seeking a non-trivial solution in the form 



such that 

W<1 

for stability. 

Substituting (9.52) into (9.51) and eliminating y, one obtains 



(9.52) 


(9.53) 


(9.54) 


where 


<5(1 + <5(co c )U + i) 2 4>i(<5 2 {co c co 2 — 8rj(o 2 ) + 2(a) c ( 1 — X) 
8\X + l) 2 (1 + <5 2 cu 2 )U + l) 2 -4X 
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In order to test the stability requirement (9.53) on the characteristic equation, i.e., 
det|J| = 0, one transforms |A| < 1 into the entire left-hand plane of the z-plane by 

A = , | A| < 1 o Re(z) < 0. (9.55) 

Carrying out the necessary algebra we have from det|J(z)| = 0 the following z- 
polynomial equation: 

(<5 3 Cco f a> 2 - 5 2 rj(o 2 + l)z 2 + (<5 (co c )z + S 2 (rja) 2 4- co 2 ) = 0. (9.56) 

A test of the polynomial equation (9.56) for possible positive real roots by the 
Routh-Hurwitz criterion (Ganthmacher [29]) indicates that the partitioned procedure 
as applied to the model coupled equations (9.47) and (9.48) give a stable solution 
provided 


(d 3 Cco c o) 2 — S 2 r](o 2 + 1) > 0. (9.57) 

Note that, if there is no position feedback (i.e., t] = 0), the model interaction equations 
solved by the present partitioned solution procedure (9.46) yields unconditionally 
stable solutions as (9.57) is automatically satisfied. Hence, a more critical stability 
assessment can be made by assuming no velocity feedback (i.e. £ = 0) for which we 
have for stability from (9.57) 



(9.58) 


The preceding stability analysis on the model interaction equations permits us to 
make the following observations. First, equation (9.58) indicates that feedback 
frequency (co c ) and the strength of the position feedback (rj) dictate the computational 
stability and not the structural frequency (cu). In other words, the position feedback 
dictates the allowable step size for stability. Thus the highest frequency of the 
controller governs stability, not the highest frequency of the structure. Since most 
controllers are designed with reduced order structure models that ignore high 
frequency dynamics, the present solution procedure is not unduly restricted by 
stability. Second, if velocity feedback is present, the allowable step size for stability 
increases until £ > >/(4^ 3 /27), at which point the solution becomes unconditionally 
stable. 

It should be noted that, instead of the stabilized form of control force equation 
(9.46d) or (9.49b), if the scalar form of (9.48) is used in the preceding stability 
analysis, the resulting stability limit is given by 


h < min 




(9.59) 
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Assuming £ « 1, the first term in the above condition allows a sufficiently large step 
size. However, since £/fj ~ 1 for a balanced control law, it imposes a step-size 
restriction h 2/cu f , which approaches the limit imposed on by a typical explicit 
integration formula. This proves the advantage of the present stabilized partitioned 
solution equation (9.46) solely from the computational stability viewpoint. 

Although not elaborated herein, a stability analysis that includes an observer 
model and the state estimation error equation has been conducted with the following 
parameter choices: 


4 = Un hzl 



0 

1 


(9.60) 


in conjunction with the structural model and the controller model already used in 
(9.49). The analysis result yields the following step-size restriction: 


h < min 



(9.61) 


It should be noted that / 2I corresponds to the Kalman filter gain magnitude which 
can be adjusted to be sufficiently small compared with (O 2 as can be assessed from 
equation (9.39b). Hence, provided / 21 < the condition given by (9.58) is seen to 
govern the maximum stable step size by the present partitioned solution procedure. 

For the general multidimensional case governed by (9.46), one observes that the 
stiffness proportional control force in practice reaches only a fraction of the total 
internal force (u = y\ Kq, f] « 1). Hence, even for a distributed stiffness proportional 
control configuration where co c co^, the stable step size given by (9.58) should 
be much larger than the maximum stable step size of a typical explicit integration 
algorithm (say, h^ x <, 21 G)^). Therefore, the computational efficiency of the present 
partitioned solution procedure is established. 


9.4 SOLUTION METHODS FOR COUPLED THERMAL-STRUCTURAL 
ANALYSIS 

Coupled thermal-structural problems are becoming a major challenge in many 
engineering disciplines such as supersonic planes, satellites, superelectronic chips, and 
jet and combustion engines. Following the finite element formulations proposed by 
Wilson and Nickell [30], Nickell and Sackman [31], Oden [32], and Oden and 
Armstrong [33], among others, the semidiscrete coupled thermal— structural governing 
equations can be written as 

Mu 4- Du + Ku — CO - f 

Q0 + H0 + e°C T u = r (9 ' 62) 

where M, D and K are the mass, damping and stiffness matrices, and f is the 
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prescribed structural loading vector; Q, H and C are the heat capacitance, heat 
diffusion and thermal expansion coupling matrices, and r is the external heat source, 
respectively; and 8 0 is the reference temperature. 


9.4.1 Conventional Implicit solution procedures 

Suppose we are given two software modules, a structural analyser and a thermal 
conduction transient analysis module and are tasked to perform the coupled response 
analysis given by (9,62). The simplest way is then to move the coupling terms CO 
and C T u in the above equation to the right-hand sides and treat them as if they are an 
applied force and an additional source term, respectively. This will permit the use of 
two single disriplined-oriented software modules for the analysis of coupled problems. 
Computationally, this amounts to employing the following staggered solution procedure : 


Mu n+1/2 + Du” + 1/2 4- Ku" +1/2 = f" +l/2 + C0£ +I/2 

Q0"+ 1/2 + H0 rt+1/2 = r" +1/2 - 0 o C T u” +U2 


(9.63) 


where + 1/2 is the predicted temperature. It turns out that if C T u" + 1/2 is predicted instead 
of 0" + 1/2 , one ends up with the same accuracy and stability limits (Park 
et al [22]), 

While the above implicit-implicit staggered procedure is simple to implement, it can 
be shown that it is only conditionally stable, even though the implicit integrators used 
to integrate the left-hand sides of (9.62) are algorithmically unconditionally stable. The 
stabilization procedure that we will describe is a mid-point rule modification of Farhat et 
al. [34]). 


9.4.2 Stabilization of Implicit-Implicit staggered solution procedure 

Stabilization of a general staggered solution procedure for coupled-field problems can be 
accomplished either by a differential-level stabilization and algebraic-level stabilization. 
In the past both stabilization strategies have been employed for fluid-structure problems, 
coupled pore fluid-soil interactions, and structure-structure interaction problems (Park 
et al. [22]; Park [23]; Felippa and Park [3]; Park and Felippa [2]). 

In general one can stabilize the implicit-implicit procedure by modifying both or 
just one of the two field equations. A successful stabilization is the one that minimizes 
the impact of stabilization in terms of software modification and computational 
overhead. Of several stabilization strategies studied, a concurrent adaptation of both 
differential and algebraic-level augmentations was found to yield the most attractive 
staggered procedure. We now outline the stabilization process. 

First, we employ a mid-point version of the trapezoidal rule as 
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yn+l/2 = y" + Sy n + UZ 

ytl + Ml _ y« + Sy n+UZ 


= 2y n + 1/2 
= 2y" + 1/2 



(9.64) 


where y can be either the displacement or the temperature vector in (9.62) and <5 is 
one-half of the step size, 5 = 1/2 At, and At is the time step size. 

Second, time discretization of the thermal coupled equations (9.62b) to obtain 

(Q + <5H)0 fl + ,/2 = <5r” + 1/2 + Q0" - <50 o C T u" +1/2 . (9.65) 


Note that in the above difference equation, the unknown structural coupling term is 
associated with the velocity u n_M/2 . It is this vector that has been found to play a 
key role in stabilization of the present procedure. In order to stabilize the extrapolation 
of the coupling term, we utilize an integrated form of u" +I/1 from the structural 
equation (9.61a): 


i" +1/a = B[Mu" + 5({ n+l/1 - Ku"+ 1/2 + C0" +I/2 )] 
B = (M + <5D)-\ 


(9.66) 


Upon substituting the above expression into (9.65), we obtain 
G0” +1/2 = R" +1/2 + d 2 0 o C T BKu£ +I/2 

rm + 1/2 = fan* 1/2 + Q 0 n _ SQ 0 C r B(Mil H + Sf”* 1 ' 2 ) 

G = Q + <5H + <5 2 0 O C T BC 


It is observed that the solution matrix G for the thermal equation is augmented with 
the additional matrix <5 2 0 O C T BC and the prediction of u£ +1/2 is replaced by 
<5BKu" + 1/ \ Also, note that the predictor for u” +1/2 is simply the previous step 
solution which has been found the most stable predictor when used in con- 
junction with the trapezoidal rule (Park [23]). 

Once the thermal equation is stabilized as described above, the structural equation 
(9.62a) can be integrated in an existing structural analysis program as if the term CO 
is an external force at each integration step. Time discretization of (9.62a) by the 
mid-point rule (9.64) yields 


Eu rt+l/2 = F" +1/2 + <5 2 C0" +1/2 

(9.68) 


E = M + <5D + <5 2 K 
F M+1/2 = <5 2 f"+ 1/2 + M(u” + <5u") + <5Du". 
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The updating procedure for states at time step (» + 1) fa achieved as follows: 

u? +I = 2u" +I ' 2 — u" 

if, + l = 2(u n+1 ' 2 - u")/<5 - u" 

^r i = 20" + 1/2 - 9" 

u n+1 = M _I (f" +1 + C0? +1 - Du; +I - Ku" + *) 

u n+ 1 = i" + <5(ii" + 1 + u") (9-69) 

u" +I = u " + <5(u" +l + 

d" +I = Q-'(r" +I - 0 o C T ii" +I - H0; +I ) 

0" +I = 0’ + d(0" +l + 0"). 


9.4.3 An analysis of stability and accuracy of stabilized procedure 

adonHnl ^ f 3ggered P rocedure Presented in (9.67)-(9.69) can be assessed by 
£? *5* l^ y T P ro “ dure - 1 for «ampfe. outlined in Park [4]. Hr* we assume 

Aa^rterized V P ^ * unif ° ra ' step inte S ratlon ™ I* 


y" +I = Xy. 

Hence, computational stability is maintained if 

Ul < I. 


(9.70) 


(9.71) 


LTuI < iZt e -Hunvitz criterion, we map the stable 

iWomarion: 0< ‘ Z ' P,ane by lhe f °' loWin 8 ****** 


u 

li " +1 

u" +I 

0” +1 

gn+i 


1 + z 
1 — z 


u 

u" 

u" 

0 n 

0 n 


(9.72) 


Substitution of (9.70) and (9.72) into (9.67)-(9.69) with D = 0 yields 

r + — <5 2 C Ifu-l fol 

|_- (I - z 2 )5 2 0 o C T M-’K z 2 Q + zSH + 5 2 0 o C T M- , cJ[_0 B J ~ \_oj (9 ' 73) 

whose characteristic equation is obtained from 
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Det| Mz 3 + VM<5z 2 + 5 2 (K + e o CQ^C T + #0 o CQ” l C T M- I K)* 

+ (5 3 VK| = 0 ( 9 - 74 ) 

where 

V = CUC T 

U = Q^H(C T Cr\ 

The reader may find a complete stability analysis in Farhat et al [34]. Hence, we 
offer the following synopsis. First, for a two-degree-of freedom problem, we have 

M = 1, K = a> 2 , Q = q, H = h, C = c (9.75) 

which, when substituted into the above characteristic equation, gives 

fljZ 3 + a 2 z 2 + a x z + a 0 = 0 (9.76) 

where 


Qr 2 

a 3 — 1, a z = 5q, a T = S 2 [co 2 + - 2 — (1 + <5 2 co 2 )J, a 0 = <5 3 aco 2 . 

Since <5, h, q, co 2 , 9 y c 1 , and m > 0, all the coefficients of the polynomial (9.76) in z 
are positive. Moreover, the quantity 

a i a i ” = 0o^ lc2< 5 3m ^ 2 (l + <5 2 co 2 ) 

is also positive, which demonstrates that the stabilized staggered solution procedure 
is unconditionally stable for the 2-d.oi. model problem. 

For multi-dimensional cases, the limiting case of K = 0, which gives rise to a 
quadratically growing structural response due to thermal coupling, can be used as a 
pathological test: 

|Mz 2 4* <5 VMz + <5 2 0 O CQ“ I C T | = 0. (9.77) 

Since M is positive definite and the other two matrices are at least semi-definite, the 
stabilized staggered procedure for this limiting case is unconditionally stable via 
Bellman's theorem [35] as successfully utilized in Park [23]. Hence, we conclude that 
the procedure given by (9.67)-(9.69) is unconditionally stable. 

As for accuracy, it can be shown that the stabilized staggered procedure is second- 

order accurate. This can be done first by expanding the difference equations (9.67)- 
(9.69) and using the governing coupled semidiscrete equations (9.62) where needed. 
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9.4.4 Computational sequence 

When we ignore structural damping, diagonal structural mass and diagonal capacitance 
matrix, the stabilized computational procedure can be summarized as follows: 


R" +I/2 = S^ 1 ' 1 + Q0" - <50 o C T B[u” + - Ku")J 


(9.78) 


(Q + <5H + <5 2 0 O C T M- , Qd" +,/2 = R" +,/2 

P +1/2 = <5 2 f" +I/2 + M(u" + <S«i") + <5Du" + <5 2 C0" +1/2 

(M + <5 2 K)u" +I/2 = p +1/J 

u n , +1 = 2u" + 1/2 - u" 

0; +l = 20 n+I/2 - e - 

u B+I = M _1 (f ,,+1 + C0? +1 — Ku? +I ) 
a n+i =i" + ^(u n+i + g") 
u B+I = u" + «5(u" +1 4- u") 
o n+l = Q~V +1 - 0 o C T u" +I - H0- + «) 

0n+I = Qn + ^- +1 + ^ 


(9.79) 

(9.80) 

(9.81) 


(9.82) 


The key for the efficiency of the above procedure compared with other candidate 
procedures is to utilize the matrix C T M -I C, which appears in the left-hand side of 
(9 ' 79 _\ is a symmetric banded matrix. Other possible stabilization involves 
CQ *C T into the left-hand side of (9.79), which has much larger bandwidth than 
the former. 

Some two-dimensional solutions of the thermal-structural interaction problems 
based on the above procedure are reported in Farhat el al. [34], * 


9.5 APPLICATION EXAMPLES 

In the preceding sections three computational methods for performing coupled-field 
dynamics analyses have been surveyed. It is anticipated that as the analyst demands 
more realistic models, all the single-field components, viz., structures, control, thermal 
and mulhbody systems, may have to be included in a typical analysis. An example 
would be a satellite undergoing solar panel deployment as well as attitude stabilization 
and vibration control. Two examples that we include are a scenario of shuttle-based 
assembly of the space station for one-cargo segment, and a vibration control of a 
generic Earth-observing platform when it is subject to a reboosting thrust. 
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MB-1 


MB-2 



9.5. 1 Manoeuvring of the SUMS with prescribed motion constraints 

Figure 9.1 illustrates the assembly of the first and second modules of the space 
station to be deployed and assembled. Each module is lifted by the shuttle remote 
manipulator system (SRMS) from the shuttle cargo bay and deployed for the first 
module and subsequently assembled into the partially assembled space station. In 
order to simulate the assembly process, first, we have studied the effect of the SRMS 
dynamics due to the required manoeuvring constraints. This incremental in-space 
construction of the space station must meet stringent geometry, weight and stiffness 
requirements as shown in Figure 9.2. The arm boom assemblies comprise two thin- 
walled graphite-epoxy circular sections called the upper arm, lower arm and end 
effector. These arms are connected by a shoulder joint (modelled by a universal 
joint), an elbow joint (modelled by a revolute joint) and a wrist joint (modelled by 
a spherical joint). The properties of the joints and arms are shown as follows (Hunter 
et al. [36]): 


(1) Upper arm: 

Young's modulus: 

E u = 1.27 x 10 n Pa 

Shear modulus: 

G u = 3.18 x 10 ,o Pa 
Length: L u = 6.38 m 

Cross-section area: 

A, = 0.0022m 1 2 * 

Moment of intertia: 

l a = 3.16 x I0 -5 m 4 

Weight: W u = 24.97 kg 


(2) Lower arm: 

Young's modulus: 

E, = 1.09 x 10" Pa 
Shear modulus: 

G, = 3.30 x 10 10 Pa 
Length: L, = 7.06 m 
Cross-section area: 

A = 0.0015 m z 
Moments of inertia: 

I, = 2.19 x 10 -5 m 4 * * * * * 
Weight: W, = 24.06 kg 
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Properties of SRMS: 

• Weight = 410 Kg 

• Length = 15 m 

• Cross Section Area = 0.0022 m* 

• Young's Module = 1.27 X 10" Pa 

• Shear Module = 3.18 X 10 10 Pa 

• Density = 1.2 X 10* Kg/m 3 * 5 

• Tip Maneuvering Speed (without payload) = 0.6 m/s 



Figure 9.2 Shuttle remote manipulator system [34], 


(3) End effector: 

Length: L e = 1.82 m 

End effector weight: 
W c = 107.14 kg 


(4) Joint weights: 

Shoulder joint weight: 
W, = 117.13 kg 
Elbow joint weight: 
Wj = 53.12 kg 
Wrist joint weight: 

W„ = 84.44 kg. 


effect of the motion constraints to the orbital motion stability can be assessed 
by modelling (1) both the space shuttle and the SRMS modeled to be rigid, (2) the 
shuttle to be rigid and the SRMS as a flexible beam (discretized into 4 elements and 
5 nodal points). By imposing angular velocity (a cubic type) at the tip of the SRMS 
(Figure 9.3), the manipulator will slew through 90° with respect to the space shuttle. 
Figures 9.4 and 9.5 illustrate the pitching angles of the rigid and flexible SRMS, and 
Figure 9.6 shows the angular velocity of the rigid and flexible SRMS. Note that the 
terminal velocities of the flexible case are non-zero, implying that the SRMS 
manoeuvring would trigger vibrations on the space shuttle modules after assembly. 
To overcome this difficulty, a more refined SRMS manoeuvring motion is necessary 
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Time, (h * 0.0001 tec) 

Figure 9.6 Response of yawing motion for flexible and rigid cases. 
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Time, see 


Figure 9,7 Imposed translational and angular motions at the tip of SRMS. 

as described in Chiou et al. [37], By adopting the refined starting and stopping 
conditions, the in-space construction of the space station can be divided into the 
following stages: 

1. manoeuvring of the SRMS to the position where its end effector is ready to 
attach the space structure which is lying in the shuttle cargo bay; 

2. contact/impact when the end effector of the SRMS collides with the space 
structure; 

3. manoeuvring of the SRMS with the space structure to attach to another space 
structure which is floating in space; 

4. contact/impact when the SRMS with the space structure collides with another 
space structure in space. 

For the first stage, the motion constraints for the tip of the SRMS are given by 
Figure 9.7 where 25 seconds of manoeuvring time is used to place the end effector 
of the SRMS to the position where the space structure/payload is located. As 
indicated in Figures 9.8 and 9.9, the angular velocity vectors for the upper arm and 
— l° wer arm of the rigid and flexible SRMS experience almost the same behaviour in 
terms of trends and magnitudes which prove that the present motion constraints are 
valid in manoeuvring the rigid and flexible SRMS. At the second stage, where the 
contact/impact has occurred, the end effector of the SRMS is approaching the 




Velocity, rad/a Velocity, rad/i 



Time* see 

Figure 9.8 Angular motions ot lower 



Time, sec 

Figure 9.9 Angular motions of upper arms. 








9.5 APPLICATION EXAMPLES 


289 


3 Contact Velocity of the End Effector and the Structure 



25 25.2 25.4 25.6 25.8 26 26.2 26.4 26.6 26.8 27 

Tune, sec 

Figure 9.10 Contact velocity and acceleration of the end effector. 


structure with velocity equal to -0.01 m/s (Figure 9.10(a)). When two bodies make 
contact at 25 s, the velocity of the end effector drops from —0.01 m/s to 0.0018 m/s 
to almost Om/s in less than one second of contact/impact time. From Figure 9.10(b), 
the contact/impact provides a peak acceleration ( — 2.4 m/s 2 ) on the end effector 
which eventually dies down because of the large mass ratio between the end effector 
and the structure. At the third stage, the SRMS lifts the structure with a motion 
constraint that is given by Figure 9.11. The purpose of this motion constraint is to 
manoeuvre the structure into the position where the previously existing structure is 
located so that the assembly of two structures can take place via contact/impact. 

From Figures 9.12 and 9.13, even though the angular velocities of the flexible SRMS 
still maintain the trends as in the rigid SRMS case, the high vibration modes can 
easily be seen as the stopping conditions of the flexible SRMS are applied. 
Consequently, due to these vibrations, the non-zero terminal velocities have occurred, 
which makes the assembly of the two structures very difficult to carry out. In 
conclusion, since the current motion constraints cannot provide the zero terminal 
velocities for the flexible SRMS, the control strategy in damping out these vibrations 
needs to be studied in order to proceed to the final stage of the present construction 
process. In Figure 9.14, the contact/impact of the fourth stage has been carried out 

by using the rigid SRMS model that the velocity of the approaching structure is 

—0.01 m/s which produces of accelerations for both assembling structures during 
1.8 s of contact/impact time. Note that after two seconds of contact/impact, both 
structures are traveling with the same velocity as indicated in Figure 9.14(a). 
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Figure 9.13 Response of upper arm during third stage. 
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Figure 9.14 Contact velocity and acceleration during third stage. 
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9.5.2 Control of earth-observing platform 

The computational efficiency achieved by the partitioned solution procedure as 
compared with the conventional solution method is sketched out in Figure 9.15. 
Assumed in the construction of the chart are the ratio of the bandwidth and the 
stiffness matrix to be 0.2, the number of the actuator and that of the sensor to be 
the same, and the ratio of the structural degrees of freedom and those of the actuator 
to be 0.1. Note that, given a real-time processing computer that can perform a wall- 
clock rate of 200 samples/second command and control, the conventional method 
can at most handle the real time control of a simple beam articulation, whereas the 
partitioned method can handle the real-time control of complex truss-beam vibrations, 
or non inear problems, the advantages of the partitioned method is more pronounced, 
as can be seen from the chart. 

The partitioned CSI simulation procedure as derived in (9.46) has been 
implemented as a stand-alone package (Park et al. [381). The present software 
implementation emphasizes the use of the widely available sequential and paralled 
analysis modules specially developed for the solution of structural dynamics equations 
Note that the solution algorithm for both the structural system and the state estimator 
is the same hence the software module, provided the right hand terms are treated 
as applied forces. Although the stabilized form of the controller and the filtered 
measurements are solved in a coupled manner, their size in general is substantially 
smaller, typically a fraction of the size of the structural system for large-scale 
problems. 

Figure 9.16 illustrates a test-bed evolutionary model of an Earth-pointing satellite. 
Eighteen actuators and 18 sensors are applied to the system (see Figure 9.16 for 
their locations) for vibration control and their locations are provided in Tables 9 1 
and 9.2 Figures 9.17-9.19 are representative of the responses for open-loop, direct 
output feedback, and dynamically compensated case does drift away initially even 
though the settling time is about the same as that by the direct output feedback 
case. However, the sensor outputs are assumed to be noise-free in these two numerical 
experiments. Further simulations with the present procedure should shed light on 
the performance of dynamically compensated feedback systems or large-scale systems 
as they are computationally more feasible than heretofore possible. 

Tables 9-3 and 9.4 illustrate the computational overhead associated with the direct 
output feedbadc vs. the use of a dynamic compensation scheme by the output present 
Kalman filtering equations, compared in those tables are for two simpler tests cases, 
viz a 3-d.oz. system and a truss beam model. In the numerical experiments herein, 
we have reiiec on the Matlab software package for the synthesis of both the control 
law gains araz the discrete Kalman filter gain matrices. Results of the full state 
feedback (FSFB) utilizing a direct feedback vs a dynamically compensated feedback 
based on the Kamian filter (K. Filter) indicate that they become competitive as the 
size of the model increases. Reported in those two tables are also the effects of 
various implementation versions, from a nominal code (version Al) to a fully 
parallelized version on a shared memory machine (allied with 8 processors). The’ 
present numerical results indicate that CPU requirements for dynamically compensated 
CSI Simula tiers would in general require about three times that of a typical structure- 
only transiemr analysis. It should be mentioned that even though there is a slight 





294 


TIME INTEGRATION METHODS FOR SYSTEM DYNAMICS 


EARTH POINTING SATELLITE DESIGN PROBLEM 



Tabic 9.1 Actuator placement for EPS 
example problem. 


Actuator 

Node 

Component 

1 

97 

X 

2 

97 

z 

3 

96 

X 

4 

96 

z 

5 

65 

y 

6 

68 

y 

7 

59 

y 

8 

62 

y 

9 

45 

y 

10 

45 

z 

11 

70 

y 

12 

70 

z 

13 

95 

X 

14 

95 

y 

15 

95 

z 

16 

95 

4>, 

17 

95 

A 

18 

95 

A 


increase of total CPU units from the compiler optimized sequential run to the parallel 
case, the actual run-clock time on the parallel machine is about one-eighth of the 
sequential case. 

To conclude, it is seen that the use of the present second-order discrete Kalman 


Deflection 
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Table 9.2 Sensor placement for EPS example 
problem. 


Sensor 

Type 

Node 

Component 

1 

Rate 

97 

X 

2 

Rate 

97 

z 

3 

Rate 

96 

X 

4 

Rate 

96 

z 

5 

Rate 

65 

y 

6 

Rate 

68 

y 

7 

Rate 

59 

y 

8 

Rate 

62 

y 

9 

Rate 

45 

y 

10 

Rate 

45 

z 

11 

Rate 

70 

y 

12 

Rate 

70 

z 

13 

Position 

95 

X 

14 

Position 

95 

y 

15 

Position 

95 

z 

16 

Position 

95 

4> z 

17 

Position 

95 

4> v 

18 

Position 

95 



EPS7 Model: Open Loop Transient Response 



— Time, _ sec 

Figure 9.17 Open loop Transient response. 
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Table 9.3 CPU results for versions of ACSIS. 


Model 

Problem 

type 

(Al) 

Nominal 

code 

(A2) 

Compiler 

optimized 

(A3) 

Parallel 

observer 

3DOF 

Transient 

6.6 

2.1 

2.1 

Spring 

FSFB 

8.0 

3.3 

3.3 


K. Filter 

12.3 

3.5 

3.3 

54 DOF 

Transient 

78.2 

5.7 

5.6 

Truss 

FSFB 

97.1 

9.4 

10.2 


K. Filter 

170.7 

13.0 

10.7 

582 DOF 

Transient 

3506 

98.6 

100.3 

EPS7 

FSFB 

7040 

190.2 

294.5 


K. Filter 

n/a 

284.2 

312.5 


Table 9.4 CPU results for ACSIS with EBE computations. 


Model 

Problem 

type 

(A4) 

E-B-E 

computation 

(A5) 

Parallel 

E-B-E 

(A6) 
Parallel 
Obs. k EBE 

3 DOF 

Transient 

3.8 

3.3 

3.3 

Spring 

FSFB 

4.9 

4.4 

4.9 


K. Filter 

6.6 

5.6 

5.0 

54 DOF 

Transient 

31.7 

13.0 

13.0 

Truss 

FSFB 

35.5 

16.9 

35.6 


K. Filter 

62.6 

27.3 

36.2 

582 DOF 

Transient 

391.7 

153.9 

n/a 

EPS7 

FSFB 

485.9 

245.9 

n/a 


K. Filter 

n/a 

n/a 

n/a 


filtering equations for constructing dynamically compensated control laws add 
computational overhead, but is only the equivalent of open-loop transient analysis 
of symmetric sparse systems of order N instead of 2N x 2N dense systems. 


9.6 CLOSING REMARKS 

The present survey have focused mostly on the research activities on the comput- 
ational methods for multibody dynamics, control-structure interactions, and coupled 
thermal-structural transients undertaken by the researchers at the Center for Space 
Structures and Controls, University of Colorado. No attempt has thus been made to 
include many important advances made by researchers around the world. A notable 
omission in this survey is computational methods for parallel computations, which 
is intensely pursued by many researchers. We hope to report on the progress on 
parallel computational methods at a later occasion. 
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Abstract 


A second-order form of discrete Kalman filtering equations is proposed as a candidate state estimator for 
efficient simulations of control-structure interactions in coupled physical coordinate configurations as opposed 
to decoupled modal coordinates. The resulting matrix equation of the present state estimator consists of the 
same symmetric, sparse N x N coupled matrices of the governing structural dynamics equations as opposed 
to unsymmetric 2N x 2 N state space-based estimators. Thus, in addition to substantial computational 
efficiency improvement, the present estimator can be applied to control-structure design optimization for 
which the physical coordinates associated with the mass, damping and stiffness matrices of the structure are 
needed instead of modal coordinates. 
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Introduction 


Current practice in the design, modeling and analysis of flexible large space structures is by and large 
based on the finite element method and the associated software. The resulting discrete equations of motion 
for structures, both in terms of physical coordinates and of modal coordinates, are expressed in a second- 
order form. As a result, the structural engineering community has been investing a considerable amount of 
research and development resources to develop computer-oriented discrete modeling tools, analysis methods 
and interface capabilities with design synthesis procedures; all of these exploiting the characteristics of 
second-order models. Recent work in the area of structural dynamics simulation and massively-parallel 
processing also rely on the second-order equation forms. 

On the other hand, modern linear control theory has its roots firmly in a first-order form of the governing 
differential equations, e.g., Kwakernaak and Sivan 1 . Thus, several investigators have addressed the issues of 
interfacing second-order structural systems and control theory based on the first-order form . As a result 
of these studies, it has become straightforward for one to synthesize direct state feedback based control laws 
within the framework of a first-order control theory and then to recast the resulting control laws in terms of 
the second-order structural systems. 

Unfortunately, controllers based on a first-order state estimator are difficult to express in a pure second-order 
form because the first-order estimator implicitly incorporates an additional filter equation 7 . However a recent 
work by Juang and Maghami 8 has enabled the first-order filter gain matrices to be synthesized using only 
second-order equations. To complement the second-order gain synthesis, the objective of the present paper 
is to develop a second-order based simulation procedure for first-order estimators. The particular class of 
first-order dynamic compensation chosen for study are the Kalman Filter based state estimators as applied to 
second-order structural systems. The proposed procedure permits simulation of first-order estimators with 
nearly the same solution procedure used for treating the structural dynamics equation. Hence, the reduced 
size of system matrices and the computational techniques that are tailored to sparse second-order structural 
systems may be employed. As will be shown, the proposed procedure hinges on discrete time integration 
formulas to effectively reduce the continuous time Kalman Filter to a set of second-order difference equations. 

The primary goal of the proposed procedure is the incorporation of this general form of state estimation 
as a simulation tool in partitioned control-structure interaction (CSI) analyses. It is expected that Kalman 
Filters for real-time control of linear time-invariant systems would be implemented in the most efficient form 
available, typically a real mode-decoupled state space realization. For analytical studies of CSI systems, 
however, where the objective is frequently simultaneous optimization of controls and structures as in Belvin 9 , 
the use of such a modal form must be weighed against the preprocessing tasks required to generate the model. 
In these cases, much more flexible controllers expressed in terms of the physical coordinates instead of the 
modal coordinates are sought, ones which can readily adapt to iterative changes in the structural parameters. 
One such control law synthesis has been proposed and demonstrated to be effective for CSI optimization . 
These studies did not account, however, for dynamic compensation when full state feedback control was 
utilized. With the discrete Kalman Filter proposed herein, a general form of dynamic compensation can be 
integrated into CSI simulation and optimization which does not impose limits on the designs of the feedback 
gains or the filter gains. 

The paper first reviews of the conventional first-order representation of the continuous second-order structural 
equations of motion, in which the state variables are defined as the displacements variables x of the second- 
order structural model and the velocities i. An examination of the corresponding first-order Kalman filtering 
equations indicates that, due to the difference in the derivative of the estimated displacement ( jp*) and the 


2 



estimated velocity (i), transformation of the first-order estimator into an equivalent second-order estimator 
requires the time derivative of measurement data, a process not recommended for practical implementation. 

Next, a transformation via a generalized momentum is introduced to recast the structural equations of motion 
in a general first-order setting. It is shown that discrete time numerical integration followed by reduction 
of the resulting difference equations circumvents the need for the time derivative of measurements to solve 
Kalman filtering equations in a second-order framework. Hence, the Kalman filter equations can be solved 
using a second-order solution software package. 

Subsequently, computer implementation aspects of the present second-order estimator are presented. Several 
computational paths are discussed in the context of discrete and continuous time simulation. For continuous 
time control simulation, an equation augmentation is introduced to exploit the symmetry and sparsity of 
the attendant matrices by maintaining state dependent control and observer terms on the right-hand-side 
(RHS) of the filter equations. In addition, the computational efficiency of the present second-order filter as 
compared to the first-order form is presented. 


Continuous Formulation of State 
Estimators for Structural Systems 

Linear, second-order discrete structural models can be expressed as 

Mi + Dx + Kx — Bu + Gw , x(0) = xo , x(0) = io (1) 

u = —Z\x — Zjx 


with the associated measurements 


Z = HyX + H 2 i + v 00 

where Af, D t K are the mass damping and stiffness matrices of size (N x N)] x is the structural displacement 
vector, (N x 1); u is the active control force (m x 1); B is a constant force distribution matrix (N x m); z 
is a set of measurements (r x 1); H x and H 2 are the measurement distribution matrices (r x N ); Z x and 
Z 2 axe the control feedback gain matrices (m x N)\ w and v are zero-mean, white Gaussian processes with 
their respective covariances Q and R\ and the superscript dot designates time differentiation. In the present 
study, we will restrict ourselves to the case wherein Q and R are uncorrelated with each other and the initial 
conditions xo and xq are also themselves jointly Gaussian with known means and covariances. 

The conventional representation of (1) in a first-order form is facilitated by 


{ x x -x 
x 2 = X = x x 

Mi 2 = Mi = Bu + Gw — Dx 2 — Kx x 

which, when cast in a first-order form, can be expressed as 

f Eq = Fq + Bu -f Gw, q — ( x x x 2 ) T 
\ z = Hq + v 


where 



0 

M 



( 3 ) 


( 4 ) 
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M 

:}■ m;> 

(5) 

It is well-known that the Kalman filtering equations 1011 for (4) can be shown to be (see, 

e.g., Arnold and 

Laub 3 ): 

Eq = Fq 4 Bu 4 EPH T R’ l z 

(6) 

where _ _ _ _ 

z — z — Hq , P — 

» ij- ’414} 

(7) 


in which U and L are positive definite matrices, q is the state estimation vector, and the matrix P is 
determined by the Riccati equation 1,3 


EPE t = FPE t + EPF t - EPH t RT 1 HPE t 4 GQG t (8) 

The inherent difficulty of reducing the first-order Kalman filtering equations given by (6) to second order 
form can be appreciated if one attempts to write (6) in a form introduced in (3): 

{ a) i\ — x 

6) £ 2 = x = xi - L x z (9) 

c) Mx 7 — —D &2 — Kx\ + Bu 4- AfL 2 z 

where 

L x = (H X U 4- HiSfRT 1 , L 7 = (H X S T 4 H 7 L) T R~ l 

Note from (9b) that £ 2 P x x . In other words, the time derivative of the estimated displacement (£) is not 
the same as the estimated velocity (i); hence, £1 and £ 2 must be treated as two independent variables, an 
important observation somehow overlooked in Hashemipour and Laub 12 . 

Of course, although not practical, one can eliminate £ 2 from (9). Assuming £1 and £ 2 are differentiable, 
differentiate (9b) and multiply both sides by M to obtain 

M x 1 ~ M£ 2 4 M L x z ( 10) 

Substituting Af£ 2 from (9c) and £ 2 from (9b) in (10) yields 

Mx 1 = —D(x 1 — L\z) — Kx\ 4 Bu 4 ML 7 z 4 ML x z (11) 

which, upon rearrangements, becomes 

Mi\ 4 Di\ 4 Kx\ — Bu 4 ML 7 z 4 ML x z 4 DL x z (12) 

There are two difficulties with the above second-order estimator. First, the numerical solution of (12) 
involves the computation of £1 when rate measurements are made. The accuracy of this computation is in 
general very susceptible to errors caused in numerical differentiation of x x . Second, and most important, the 
numerical evaluation of z that is required in (12) assumes that the derivative of measurement information is 
available which should be avoided in practice. We now present a computational procedure that circumvents 
the need for computing measurement derivatives and that enables one to construct estimators based on the 
second-order model form. 
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Second-Order Transformation of 
Continuous Kalman Filtering Equations 

This section presents a transformation of the continuous time first-order Kalman filter to a discrete time set of 
second-order difference equations for digital implementation. The procedure avoids the need for measurement 
derivative information. In addition, the sparsity and symmetry of the original mass, damping and stiffness 
matrices can be maintained. Prior to describing the numerical integration procedure, a transformation based 
on generalized momenta is presented which is later used to improve computational efficiency of the equation 
solution. 


Generalized Momenta 

Instead of the conventional transformation (3) of the second-order structural system (1) into a first-order 
form, let us consider the following generalized momenta (see, e.g., Jensen 13 and Felippa and Park ). 

f a) x x =x ( 13 ) 

\ 6) X 2 — AM x\ + Cx\ 

where A and C are constant matrices to be chosen. Note that AM should be nonsingular in order to obtain 
an equivalent form of (1). Time differentiation of (13b) yields 

ir 2 = AMx\ + Cx\ (14) 


Substituting (1) via (13a) into (14), one obtains 

x 2 — A(Bu -4- Gw) — (AD — C)x i — AKxy 
Finally, pairing of (13b) and (15) gives the following first-order form: 


f AM °l/M + [ c - 7 1M = 

[AD-C I J \± 2 J + [AK o J \ * 2 j 
0 


r 0 1 
A(Bu + Gw) J 


[A(Bu + Gw)\ 

The associated Kalman filtering equation can be shown to be of the following form: 


AM 

AD-C 


oiNii.r c -miiif o \ 

/ Jl £ 2j [AK Oj\i 2 J \ ABu ) 


AM 

°1 f 

AD-C 

/][ 


b 

L 2 


where 


To- 1 


h = (H X V + H 2 S) t R-\ U = ( H\S T + H 2 L) 1 R 

and Hi and H 2 correspond to a modified form of measurements expressed as 

z = H\x + H 2 x = Hixi + H 2 x 2 


(15) 


(16) 


(17) 


( 18 ) 


where 

Hi = Hr - H 2 M- l A~ 1 C, H 2 = H 2 M~ l A~ l 
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Clearly, as in the conventional first-order form (9), x\ and x 2 in (17) are now two independent variables. 
Specifically, the case ofA = A/“ l andC = 0 corresponds to (3) with x 2 = i\. However, as we shall see 
below, the Kalman filtering equations based on the generalized momenta (13) offer several computational 
advantages over (3). 


Numerical Integration 

At this juncture it is noted that in the previous section one first performs the elimination of f 1 in order 
to obtain a second-order equation, then performs the numerical solution of the resulting equation. This 
approach has the disadvantage of having to deal with the time derivative of measurement data. To avoid 
this, we will first integrate numerically the associated Kalman filtering equation (17). 

The direct time integration formula we propose to employ is a mid-point version of the trapezoidal rule. 


a ) 

b) 



( f., i n r i, 

, .+!/! 

{* } +6 {l 

} 


M n 

U) -( 

•I} 


(19) 


where the superscript n denotes the discrete time interval t n — n/i, h is the time increment and 6 h( 2. 

It should be noted that we have chosen the trapezoidal rule due to its unconditional stability and high 
accuracy while it does not introduce any numerical damping (see, for example, Dahlquist and Park ). 
Contamination of damping from numerical dissipation can not only adversely affect the solution accuracy 
but lead to misinterpretation of the simulation results. 


Time discretization of (17) by (19a) at the n + 1/2 time step yields 


AM 

AD-C 


01 /x? +1/2 -*? 

/] u +1/2 -xs 




= 6 


AM 

AD-C 



z n+l/2 + 6 | ^ flu „ + i /2 } 


( 20 ) 


The above difference equations require the solution of matrix equations of 2 N variables, namely, in terms 
of the two variables Xj* 1 ^ 2 and x”* 1 ^ 2 , each with a size of N . To reduce the above coupled equations of 
order 2 N into the corresponding ones of order N, we proceed in the following way by exploiting the nature 
of parametric matrices of A and C as introduced in (13). To this end, we write out (20) as two coupled 
difference equations as follows: 

,4M(x" +1/2 - *?) + 6{Cx n l +in - x" +1/2 ) 

= 6AML 1 z n + 1 ' 2 (21) 

(AD - C)(x n + 1 ' 2 - x") + (x^ +1/2 - x^) + SAKx ? +1/2 

= 6(AD - C)Liz n+1/2 + 6L 7 z n+l/ 2 + 6ABu n+l/2 (22) 

Multiplying (22) by 6 and adding the resulting equation to (21) yields 

A(M + 6D + S 2 I< )x” +1/2 = (AM + S(AD - C))x" + Sx ” 

+{6AMLi + 6 2 (AD - C)L l + 6 2 L 2 }f n+l/2 + 6 2 ABu n+l/2 (23) 
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Of several possible choices for matrices A and B, we will examine the following two specific cases: 


f a) A = I , C= D 
{ b) A — M~ l , C- 0 


( 24 ) 


where the mass matrix M is nonsingular due to its physically positive definite nature since the kinetic energy 
of structural system is positive for any admissable motion. It is noted that the above two choices, although 
mathematically equivalent, lead to different computational implementations as discussed below. 

The choice of (24a) reduces (23) to: 


(M + SD + 6 2 K)xl +1/2 = Mi" + Sx’i + 6 2 Bu n+lf 2 
+6{ML t + SL 2 }z n+l/2 

so that once i" +l ^ 2 is computed, i, + l ^ 2 is obtained from (22) rewritten as 

ij +1/2 = i? + 6g n - 6Kx" +l/2 


(25) 


(26) 


where 

g n = Bu n + 1 * 2 + Z 2 r n+l/2 (27) 

which is already computed in order to construct the right-hand side of (25). Hence, Kx" +1 ^ 2 is the only 
additional computation needed to obtain x!^ 1 ^ 2 - It is noted that neither any numerical differentiation nor 
matrix inversion is required in computing x 1 ^ 1 ^ 2 . This has been achieved through the introduction of the 
general transformation (13) and the particular choice of the parameter matrices given by (24a). 

On the other hand, if one chooses the conventional representation (24b), the solution of x t + ^ is obtained 
from (23) 

(M + 6D + 6 2 K)x”+ lf2 = (M + 6D)x” + 6Mx 5 

+6{{M + 6D)h + <5ML 2 }z n+1/2 + 6 7 Bu n + l/7 (28) 

Once x^ 1/2 is obtained, x 2 + 1/f2 can be computed either by 


.n+l/2 

x 2 


= (x” +1/2 -xJ)/6-Z 1 f n+l/2 


(29) 


which is not accurate due to the numerical differentiation to obtain x t + ^ , or by (22) 


x" +1/2 = x” + Sg n - 6M~ l Kx"+ l/7 — 
M' l D(x" f 1/2 - x?) -f 6M ' 1 Z?L 1 r n+1/2 


(30) 


which involves two additional matrix-vector multiplications, when D ^ 0, as compared with the choice of 
A = I and C = D. Thus (24a) is the preferred representation in a first-order form of the second-order 
structural dynamics equations (1) and is used in the remainder of this work. 


Decoupling Of Difference Equations 

We have seen in the previous section, instead of solving the first-order Kalman filtering equations of 2n 
variables for the structural dynamics systems (1), the solution of the implicit time-discrete estimator equation 
(25) of n variables can potentially offer a substantial computational saving by exploiting the reduced size 
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and sparsity of M,D and K. This assumes that z n+1/2 and u n+1/2 are available, which is not the case since 
at the n th time step 


u n+l f 2 = -Z l x” +l,, -Z 2 fj +1/2 

(31) 

£"+!/ 2 = Z n+l/ 2 - Hxi " +1/2 - H 2 Xt +l/2 

(32) 


requires both x" +1/2 and Xj +l/2 even if z n+1/2 is assumed to be known from measurements or by solution 
of (1). Note in (32), the control gain matrices are transformed by 

Z x - Z x - Z 2 M~ l A~ l C, Z 2 = Z 2 M~ l A~ 1 


There are two distinct approaches to decouple (25) and (26) as described in the following sections. 


Discrete Time Update 

For systems utilizing discrete-time (i.e. sample and hold) control, equations (31) and (32) become 

u n+i/2~_Z 1 f^-7 2 x5 (33) 

+ ~ Z n - ( 34 ) 

The time integration step size of the estimator must then be equal to the sample rate of the control, while the 
continuous structural equations may also be integrated at the same rate or at some fraction of the sampling 
rate for simulation accuracy considerations. For the present purposes, we have assumed that the sampling 
interval is the same as the integration time stepsize. 

Discrete time simulation is quite simple to implement as the control force and state corrections are treated 
with no approximation on the right-hand-side (RHS) of (25) and (26). Should continuous time simulation 
be required, a different approach is necessary. 


Continuous Time Update 


To simulate the system given in (25) and (26) in continuous time, strictly speaking, one must rearrange 
(25) and (26) so that the terms involving x" +1/2 and x? + ' /7 are augmented to the left-hand-side (LHS) of 
the equations. However, this augmentation into the solution matrix ( M + 6D + 6 2 K) would destroy the 
computational advantages of the matrix sparsity and symmetry. Thus, a partitioned solution procedure has 
been developed for continuous time simulation as described in Park and Belvin 17 . The procedure, briefly 
outlined herein, maintains the control force and state correction on the RHS of the equations as follows. 

First, x? +1/2 and i 2 +l ^ 7 are predicted by 


-n + l/2 
Z 1 p 


'1 > 


-n + l/2 
x 2p 


(35) 


However, instead of direct substitution of the above predicted quantity to obtain «p + 1/2 and z£ +1/2 based 

. , , n + l/2 . -n + l/2 

on (31) and (32), equation augmentations are introduced to improve the accuracy ot u p ana z p 

Of several augmentation procedures that are applicable to construct discrete filters for the computations of 
u n +i/2 an< j .j n +i/2, we substitute (26) into (31) and (32) to obtain 


' „"+‘/2 = -Zxi" + l/7 - Z 2 {*2 - 61<i n x + l,2 + 

6Bu n+1 / 2 + SL 2 z n+l f 2 ) 

* f" + 1 / 2 = 2" + 1/2 - ^li7 +1/2 - 

f/ 2 (x? - SKx ? +1/2 + 6Bu n+l f 2 + 6L 2 z n+1 / 2 ) 
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Rearranging the above coupled equations, one obtains 

\(I + SZ 2 B) 6Z 2 L 2 1 f u n+1/2 

[ 6H 2 B (I + 6H 2 l 2 )\ \f n+1/2 

r -Z 2 5Z-(Z l -SZ t K)zl*‘ l/i \ (37) 

{ 2 «+i/2 _ H 2 i- - (ff, _ *i? 2 J<)i? +l/2 J 

which corresponds to a first order filter to reduce the errors in computing x 2 = Mi + Di. A second-order 
discrete filter for computing u and z can be obtained by differentiating u and z to obtain 

/ « = — — ^ 2*2 , (38) 

\ z — z — H ii\ — H 2 x 2 



and then substituting x x and x 2 from (17). Subsequently, (19) is applied to integrate the equations for u 
and z which yields 


{i + SZiB + S'ZtM-'B _ 

[ 6{H 2 B + 6H\M~ 1 B) I + 6Hi(Li + SM~ 1 L 2 ) + 6H 2 L 2 \ \ z n+1/2 J 

fu"l . ( Z 1 M-\xl-6Kx n l +l/7 -Dx n l +lf2 ) + Z 2 Kx n 1 +1,2 \ f 0 1 

i i" J 6 1 H x M- l {x 5 - *Kf? +1/2 Dx? +1/2 ) + H 7 Kx n l +l/2 } + l *" +1/2 J 


u{; 


(39) 


The net effects of this augmentation are to filter out the errors committed in estimating both xi and x 2 . 
Solution of (39) for tt n+1 / 2 and F l+1 ^ 2 permits (25) and (26) to be solved in continuous time for x” +1/2 and 
^n+i /2 Subsequently, (19b) is used for x" +l and xj" 1 * 1 . 

The preceding augmentation (39) leads to an accurate estimate of the control force and state estimation error 
correction at the (n+1/2) time step. Although (39) involves the solution of an additional algebraic equation, 
the equation size is relatively small ( size = number of actuators (m) plus the number of measurements (r) 
). Thus, (39) is an efficient method for continuous time simulation of the Kalman filter equations provided 
the size of (39) is significantly lower than the first order form of (4). The next section discusses the relative 
efficiency of the present method and the conventional first order solution. More details on the equation 
augmentation procedure (39) may be found in Park and Belvin 17 . 

Finally, it is noted that by following a similar time discretization procedure adopted for computing 
and Xj +1 ^ 2 , the structural dynamics equation (1) can be solved by 


/ (M + 6D + 6 2 K)x n + l/2 = Mx J 4- + 6 2 Bu n + 1/2 (40) 

\ x^ l/2 = xS + <5£u n+i / 2 - 6Kz1+ l/2 

Thus, numerical solutions of the structural dynamics equation (1) and the filter equation (20) can be carried 
out within the second-order solution context, thus realizing substantial computational simplicity compared 
with the solution of first-order systems of equations (4) and the corresponding first-order estimator equations 
( 6 ). 

It is emphasized that the solutions of both the structural displacement x and the reconstructed displacement 
x employ the same solution matrix, (Af -f SD 4- 8 2 K). The computational stability of the present procedure 
can be examined as investigated in Park 18 and Park and Felippa 19,20 . The result, when applied to the present 
case, can be stated as 

6 2 A m ax < 1 ( 41 ) 
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where A m ax is the maximum eigenvalue of 


(A 2 / + A Z 2 B + Z l M- l B)y = 0 (42) 

Typically the control laws are formulated in terms of low-frequency response components, viz., 

BotG T KG (43) 

for the displacement feedback case where G is a projection matrix that extracts only low-frequency compo- 

nents from the structural stiffness matrix. Hence, A max is in practice several orders of magnitude smaller 
than /imax of the structural dynamics eigenvalue problem: 

pMy = Ky (44) 

Considering that a typical explicit algorithm has its stability limit ax * ^ ^ the maximum step size 
allowed by (42) is in fact several orders of magnitude larger than allowed by any explicit algorithm. 


Computational Efficiency 

Solution of the Kalman filtering equations in second-order form is prompted by the potential gain in com- 
putational efficiency due to the beneficial nature of matrix sparsity and symmetry in the solution matrix 
of the second-order estimator equations. There is an overhead to be paid for the present second-order pro- 
cedure, that is, the additional computations introduced to minimize the control force and state estimation 
error terms on the right-hand-side of the resulting discrete equations. The following paragraphs show the 
second-order solution is most advantageous for estimator models with sparse coefficient matrices M, D and 
K . 

Solution of the first order Kalman filter equation (6) or the second-order form (25-26, 39) may be performed 
using a time discretization as given by (19). For linear time invariant (LTI) systems, the solution matrix is 
decomposed once and subsequently upper and lower triangular system solutions are performed to compute 
the estimator state at each time step. Thus, the computations required at each time step result from 
calculation of the RHS and subsequent triangular system solutions. For the results that follow, the number 
of floating point operations are estimated for LTI systems of order 0(N). In addition, it is assumed that the 
mass, damping and stiffness matrices (M,£> and K) are symmetric and banded with bandwidth aN i where 
0 < a < (0.5- 2^). 

The first-order Kalman filter equation (6) requires (AN 2 + 2Nr + O(N)) operations at each time step. The 
discrete time second-order Kalman filter solution (25-26, 33-34) require (8a 2 N 2 +2aN 2 +ZNm+ANr+0(N)) 
operations and the continuous time second-order Kalman filter (25-26, 39) require (8acN 2 + 2otN +5/'/m + 
6Nr + (r+ m) 2 -f O(N)) operations at each time step. To examine the relative efficiency of the first-order 
and second-order forms, several cases are presented as follows. 

First, a worst case condition is examined whereby M, D and I\ are fully populated (a = 0.5 — 2 V) and 
r = m s N . Only for this extreme condition with large numbers of sensors and actuators relative to the 
system order, the first order Kalman filter becomes somewhat more efficient than the second-order discrete 
Kalman filter solution presented herein. 

For typical structural systems, M and K are almost always banded. In addition, the number of sensors 
and actuators is usually small compared to the system order N. If the number of actuators (m) and the 
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number of measurements (r) are proportional to the bandwidth ( r = m = the second-order discrete 

Kalman filtering equations become computational attractive as long as a < 0.394. It should be noted that 
the larger the size of the structural systems, the smaller the bandwidth becomes, with the range of a to be 
0.05 < a < 0.15. 

Finally, for the special case of modal-based structural models, one has a — * 0. For this case, as long as 
sensors and actuators are sufficiently smaller than the modal degrees of freedom, the present second-order 
state estimator can be substantially more efficient than the classical first-order form. This is because the 
conventional state space-based estimator must deal with a fully coupled nonsymmetric 2N x 2 N system 
whereas the present second-order estimator deals with a diagonal N x N system. A more detailed discussion 
can be found in Belvin 9 . 


Implementation and Numerical Evaluations 

The second-order discrete Kalman filtering equation derived in (25) and (26) have been implemented along 
with the stabilized form of the controller u and the filtered measurements z in such a way the estimator 
computational module can be interfaced with the partitioned control-structure interaction simulation package 
developed previously by Belvin9, Park and Belvin 17 Alvin and Park 21 . It is emphasized that the solution 
procedure of the present second-order discrete Kalman filtering equations (25) and (26) follows exactly the 
same steps as required in the solution of symmetric, sparse structural systems. It is this attribute that makes 
the present discrete filter attractive from the simulation viewpoint. For a succinct comparison between the 
present CSI simulation procedure and conventional state space-based simulation procedures, the equations 
that need to be implemented in both of the procedures are summarized below. 

Partitioned Control-Structure Interaction Equations 

The partitioned procedure for simulating the control-structure interaction problems developed in Belvin 9 and 
Park and Belvin 17 exploits the second-order diferential equation form whenever possible as shown below. 


r Structure: a) 

Sensor Output: 6 ) 

i Estimator: c) 

Control Force: d) 

K Estimation Error: e) 


Mq + Dq + Kq = f-f Bu + Gw 

q(0) = q 0 , q(0) = q 0 

z = Hx 4 * v 


M 

0 

q( 0 ) = 0 , 


m + [“ vkih+’-h 


M 

0 


p( 0) = f(0) + Bu(0) 



(45) 


u + F 2 M~‘Bu = F 2 (M -, p + L 2 -r) + F,^ 


7 + H„L 27 =i- H„M - 1 (f> - Bu) - H d 4 


In addition, notice that the control laws (u) and the estimation error ( 7 ) are parabolically stabilized and 
solved in a separate software module from the estimator and the structural analyzers, thus effectively ren- 
dering a computaionally efficient and accurate procedure. 
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Conventional Control- Structure Interaction Equatioons 

In contrast to the partitioned procedure summarized above, conventional control-structure interaction simu- 
lation employs a first-order differential equation form as shown below, thus requiring the solution of 2n x2n)- 
system equations for structures and the observer. In addition, the control laws and the estimation errror are 
not stabilized, which can give rise to an accumulation of computational errors. 


where 


and 


r Structure: a) 

Sensor Output: b) 

Estimator: c) 

Control Force: d) 


l Estimation Error: e) 


x = Ax + Ef+Bu + Gw 
x(0) = x 0 

z — Hx -|- v 

x = Ax + Ef + Bu -F lor 
x(0) = 0 

u = — Fx 
7 = z- Hx 



(46) 


Numerical Experiments 

The first example is a truss beam shown in Fig. 1, consisting of 8 bays with nodes 1 and 2 fixed for 
cantilevered motions. Actuator and sensor locations, as well as their orientation, are given in Table 1. 
In the numerical experiments reported herein, we have relied on the Matlab software package 22 for the 
synthesis of both the control law gains and the discrete Kalman filter gain matrices. Figures 2, 3 and 4 
show the vertical displacement time response at node 9 for open-loop, full state feedback, and dynamically 
compensated feedback cases, respectively. In the present pepper, a full state feedback corresponds to the case 
for which the number of sensors are the same as the total system degrees of freedom whereas the dynamically 
compensated case corresopnds to a smaller number of sensors as compared with the total system degrees of 
freedom. Note the effectiveness of the dynamically compensated feedback case with four actuators and six 
sensors as indicated in Table 1 by the present second-order discrete Kalman filtering equations as compared 
with the full state feedback cases. 

Figure 5 illustrates a testbed model of an Earth-pointing satellite. For vibration control, 18 actuators and 
18 sensors are configured throughout the system; their locations are provided in Tables 2 and 3. Figures 6, 7, 
and 8 are a representative of the responses for open-loop, full state feedback, and dynamically compensated 
cases, respectively. In both examples, the estimator states are the estimated physical displacements and 
generalized momenta as previously developed, and thus the number of effective states is equal to 2 N, where 
N is the number of physical displacement variables of the second-order structural system. Therefore, the 
Kalman filter for the truss example has 108 states, and the filter for the satellite has 1164 states, a substantial 
increase over typical estimator orders for such systems. Further simulations with the present procedure should 
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shed light on the performance of dynamically compensated feedback systems for large-scale systems as they 
are computationally more feasible than heretofore possible. 

The computational overhead associated with the full state feedback vs. the use of a dynamic compensation 
scheme by the present Kalman filtering equations is reported in Table 4. It is seen that the use of the 
present second-order discrete Kalman filtering equations for constructing dynamically compensated control 
laws adds computational overhead, only an equivalent of open-loop transient analysis of symmetric sparse 
systems of order N instead of 2 N x 2 N dense systems. This is evidenced in Table 4 in that the normalized 
CPU time for the dynamically compensated case (designated as K. Filter) is 63.16 whereas the total CPUs 
for the full state feedback case (FSFB) plus that of the open loop dynamic response (Transient) is 64.18. 


Summary 

The present paper has addressed the advantageous features of employing the same direct time integration al- 
gorithm for solving the structural dynamics equations and for integrating the associated continuous Kalman 
filtering equations. The time discretization of the resulting Kalman filtering equations is further facilitated 
by employing a canonical first-order form via a generalized momenta. When used in conjunction with the 
previously developed stabilized form of control laws in Park and Belvin 17 , the present procedure offers a sub- 
stantial computational advantage over the simulation methods based on a first-order form when computing 
with large (i.e. nearly full system dynamics) and sparse estimator models. 

In order to minimize the deleterious effect of numerical damping and phase distortion in the solution of the 
discrete Kalman filtering equations, the trapezoidal rule is employed. This is due to the wellknown fact that 
the trapezoidal rule conserves the system energy with minimum phase error among all the time integration 
formulas of second-order accuracy 15,16 

Computational stability of the present solution method for the filter equation has been assessed based on 
the stability analysis result of partitioned solution procedures 18 . To obtain a sharper estimate of the stable 
integration step size, a more rigorous computational stability analysis is being carried out and will be reported 
in the future. 
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Table la 


Actuator Placement for Truss Example Problem 


Actuator 

Node 

Component 

1 

2 

y 

2 

18 

y 

3 

9 

y 

4 

9 

X 


Table lb 

Sensor Placement for Truss Example Problem 


Sensor 

Type 

Node 

Component 

1 

Rate 

2 

y 

2 

Rate 

18 

y 

3 

Rate 

9 

y 

4 

Rate 

9 

X 

5 

Position 

9 

y 

6 

Position 

9 

X 





Table 2 

Actuator Placement for EPS Example Problem 


Actuator 

Node 

Component 

i 

97 

X 

2 

97 

z 

3 

96 

X 

4 

96 

z 

5 

65 

V 

6 

68 

y 

7 

59 

y 

8 

62 

y 

9 

45 

y 

10 

45 

z 

11 

70 

y 

12 

70 

z 

13 

95 

X 

14 

95 

y 

15 

95 

z 

16 

95 

<f>z 

17 

95 

<f>y 

18 

95 

4>z 





Table 3 

Sensor Placement for EPS Example Problem 



3 

Rate 

96 

X 

4 

Rate 

96 

z 

5 

Rate 

65 

y 

6 

Rate 

68 

y 

7 

Rate 

59 

y 

8 

Rate 

62 

y 

9 

Rate 

45 

y 

10 

Rate 

45 

Z 

11 

Rate 

70 

y * 

12 

Rate 

70 

z 

13 

Position 

95 

X 

14 

Position 

95 

y 

15 

Position 

95 

z 

16 

Position 

95 

<t>: 

17 

Position 

95 

<j)y 

18 

Position 

95 

<h 




Table 4 

CPU Results for ACSIS Sequential and Parallel Versions 


Model 

Problem 

Type 

Sequential 

Parallel 

54 DOF 

Transient 

4.5 

5.6 

Truss 

FSFB 

9.4 

10.2 


K. Filter 

13.0 

10.7 

582 DOF 

Transient 

98.6 

100.3 

EPS7 

FSFB 

190.2 

294.5 


K. Filter 

284.2 

321.5 
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Figure Is TRUSS BEAM PROBLEM 
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Truss Model: Open Loop Transient Response 
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Figure 2: Truss Transient Response 
Truss Model: Full State Feedback Response 
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Figure 3: TRUSS FULL STATE FEEDBACK RESPONSE 







Deflection 


Truss Model: Controlled Response with Kalman Filter 
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Figure 4: Truss Response with Filter 






Figure 5: GENERIC EARTH POINTING SATELLITE EXAMPLE 
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EPS7 Model: Open Loop Transient Response 
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EPS7 Model: Full State Feedback Response 
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EPS FULL STATE FEEDBACK RESPONSE 


Model: Controlled Response w/Kalman Filter 
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Dynamics of Adaptive Structures: Design through Simulations 
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1. INTRODUCTION 

Mechanisms and static stress analyses have long been the major considerations in the 
design of many articulated structures or adaptive structures in the past. However, high- 
performance requirements on these structures have added the dynamics considerations 
as a new added design criterion in recent years. This is especially true in the design of 
adaptive or deployable space structures that involve the combined phenomena of the 
orbital mechanics, structural configuration changes and flexible vibrations in a coupled 
manner. Hence, little attention has been given to, in the design of reconfigurable flexible 
space structures, the influence of the accompanying dynamics during the maneuvering 
as an integral part of the design requirements. 

The adaptations of human bodies, animals and bacteria to spatial dynamical motions 
have been previously studied[l-3]. Recently, several investigators developed the so- 
called angular momentum preserving rotational maneuvering control algorithms and 
applied them to robotics and spacecraft attitude controls[4-6]. As a result, the intrinsic 
adaptations of the momentum conservation (violation for that matter) laws by spring 
board divers, ice skaters as well as gymnasts are well understood, which have been 
subsequently utilized for the design of space robotics maneuvering and space rendezvous 
scenarios. These studies have dealt mostly with rigid bodies linked by frictionless joints 
and focused on the development of various control algorithms for nonholonomic rigid 
dynamical system. 

The use of a helical bi-morph actuator /sensor concept [7] by mimicking the change of 
helical waveform in bacterial flagella is perhaps the first application of bacterial motions 
(living species) to longitudinal deployment of space structures. However, no dynamical 
considerations were analyzed to explain the waveform change mechamsms[3, 7]. The 
objective of the present paper is to review various deployment concepts from the dy- 
namics point of view and introduce the dynamical considerations from the outset as part 
of design considerations. Specifically, the impact of the incorporation of the combined 
static mechanisms and dynamic design considerations on the deployment performance 
during the reconfiguration stage is studied in terms of improved controllability, maneu- 
vering duration and joint singularity index. It is shown that intermediate configurations 
during articulations play an important role for improved joint mechanisms design and 
overall structural deployability. 




2. EXAMPLES OF ADAPTIVE STRUCTURES 


2.1 Bacterial Flagella 

In studying the chemotaxis of bacteria such as Salmonella , scientists discovered that 
their motions are intertwined with smooth swimming interrupted by short periods of 
tumbling[3]. In particular, the change in waveforms do not follow the intuitive way, 
vz., from one normal wave form to the adjacent discrete wave state. Instead, the 
transition of the waveform jump from one wave sometimes to its half-length wave. 
Calladine[3] conjectured that the intermittent existence of bi-stable subunits along 
the helical flagella structure are responsible for the formation of partly stable curly 
right-handed helices. It is these bi-stable subunits that cause jumps in the waveform 
formation. 




Fig. 1 Possible Waveforms of Flagella of Salmonella 

From the mechanical deployment perspectives, the large motions due to the jumps in 
waveform change in bacterial flagella pose the following questions: 1) how can such 
large motions be possible what are the sources of the torques that make such large 







motions possible?; 2) axe those motions created by minimizing the energy requirements 
or by triggering unstable motion paths so that the energy need remains minimal?; 
3) can the large motion phenomenon be explained solely by quasi-static equilibrium 
considerations or be explained only by the dynamical considerations? 

Experiments as well as analytical studies[3] so far identified twelve polymorphic helical 
forms with a tubular chains of 20 nanometer in diameter as shown in Fig. 1. 

2.2 rteconflgurable Truss Beams 

Figure 2 illustrates three representative reconfigurable truss beams. The sequentially 
deployable maneuvering tetrahedral beam is shown in Fig. 2a and can only be deployed 
sequentially, hence can’t simulate the jumps in waveform of the bacterial flagella. 



Fig. 2 Various Reconfigurable Truss Beams 

By inserting actuator-encorder pairs into some of the truss members as shown in the 
variable geometry truss (Fig. 2b), it is possible to shape the beam as desired. The 
batten actuated beam as shown in Fig. 2c is perhaps the simplest reconfigurable truss. 
In both the last two cases, the actuators may be viewed as bi-stable subunits which, 
unlike for the case of tumbling motions in flagella case, do require control forces. 






3. NONHOLONOMICALLY CONTROLLED RECONFIGURABLE STRUCTURES 


The equations 
be written as 


of motion for nonholonomicaUy controlled reconfigurable 
p = f(<) — S(q) + Bu + CA 


p = Mq + D(q) 


structures can 
( 1 ) 


with the constraints: 


*/c(q) =0 


*w(q) =° 


c = 

B = 


d<& K 

dq 

3q 


( 2 ) 


In the above equations, p is the generalized momenta, f is the applied force, S is the 
internal force, u and A are the nonholonomic and kinematic contraint forces, M is the 
generalized inertia matrix, D is the damping operator, and $*(q) and $w(q) are system 
kinematic and nonholonomic constraint equations. It should be noted that both u and 
A can be augmented with active control forces, when necessary. 


Figure 3 illustrates a design example that involves the sizing of the double moemnt 
gyros[ll] for effecting the maneuvering as well as the necessary vibration control. The 
moment gyros can in turn be made of from micro to mini sizes (12], depending upon 
the torque requirements. In this particular example, the task is to shape the artic- 
ulated straight beam to form an hexagonal polygonal structure in space or can be 
shaped to form a helix if desired. Therefore, the role of gyros is to perform triple tasks 
concurrently: maneuvering, vibration control, and if necessary bi-stable units for easy 
articulation. 





Fig. 3 Articulation of Beam-Like Structure via Control Moment Gyros 




It should be mentioned that, for three rigidly linked planar maneuvering that con- 
serves the angular momentum, the problem has been analyzed in [6]. It is for flexible 
cases the solution can be complicated. These and solutions of other related problems 
will be reported at the conference. 
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